# Run regressions for foot traffic analysis


# Setup ----
options(scipen = 999)
library(caret)
library(tidyverse)
library(lfe)
library(stargazer)
library(cragg)
library(xtable)


load("/02_Data/FL_visitor_1804_1910_tract.RData") # NOTE: Proprietary data, values of all variables have been replaced with * 
load("/02_Data/tract_flood_FL.RData")             # NOTE: Proprietary data, values of all variables have been replaced with * 
load("/02_Data/NFIP_Michael_Claim.RData")
load("/02_Data/tract2017.RData")
load("/02_Data/NFIPTakeup.RData")
load("/02_Data/FloodZone.RData")
load("/02_Data/SavingsAccount.RData")
load("/02_Data/SCI_v3.RData")


# tract-level data as of 2017 --------------------------------------------------
tract_data = tract2017 %>%
  left_join(tract_flood_FL, by=c("GEOID"="FIPS_tract")) %>%
  left_join(NFIPTakeup, by="GEOID") %>%
  left_join(FloodZone, by="GEOID") %>%
  left_join(SavingsAccount, by="GEOID") %>%
  left_join(SCI_v3, by=c("GEOID"="tract")) %>%
  mutate(FloodDepth_ft_tract_dummy = ifelse(FloodDepth_ft_tract>0, 1, 0))



# tract-business-time level data, April 2018 to October 2019 -------------------
FL_visitor_1804_1910 = FL_visitor_1804_1910_tract %>%
  left_join(NFIP_Michael_Claim, by=c("visitor_home_tract"="CensusTract", "month")) %>%
  left_join(tract_data, by=c("visitor_home_tract"="GEOID")) %>%
  mutate(
    FloodDepth_ft_tract_dummy = ifelse(FloodDepth_ft_tract>0, 1, 0),
    ClaimAmount               = ifelse(is.na(ClaimAmount), 0, ClaimAmount),
    WAvgPaidWeek              = ifelse(is.na(WAvgPaidWeek), 0, WAvgPaidWeek),
    month                     = factor(month)
  )

date.binarized = predict(dummyVars(formula = "~ month", FL_visitor_1804_1910, sep = '_'), FL_visitor_1804_1910)        # the "dummyVars" formula in the "caret" package can be used to binarize the categotical variable automatically
colnames(date.binarized) = substr(colnames(date.binarized)[1:length(colnames(date.binarized))], 7, nchar(colnames(date.binarized)[1:length(colnames(date.binarized))]))

FL_visitor_1804_1910 = cbind.data.frame(FL_visitor_1804_1910, date.binarized)
rm(date.binarized)





# ---------------------------------------
# ------ Propensity score matching ------
# ---------------------------------------
tract_data$FloodInsTakeup_High = ifelse(tract_data$FloodInsTakeup>0.5, 1, 0)

library(MatchIt)

mod_match <- matchit(FloodInsTakeup_High ~ log_medincome2017 + log_population2017 + ptg_Unemployed2017 +
                       ptg_white2017 + ptg_OwnerOccupied2017 + ptg_Mortgage2017 + ptg_Bachelor2017 + Gini2017 + 
                       ptg_SavingsAccount2017 + FloodZone_Share_Dev_Area_tract, 
                     method = "nearest", 
                     distance = "glm",
                     data = tract_data[tract_data$FloodDepth_ft_tract>0,])

tract_FL_flooded_match = match.data(mod_match)



# Matched sample --------------------------------------------------------------
FL_visitor_1804_1910_matched = FL_visitor_1804_1910 %>%                 
  filter(visitor_home_tract %in% tract_FL_flooded_match$GEOID)








################################################################################
# Table 3: Effect of Flood Insurance Take-up 
################################################################################

a0 = felm(log(1+visitor_count_home_tract) ~ (FloodDepth_ft_tract + 
                                               
                                             FloodDepth_ft_location + NAICS_2digit_name + dummy_brands +
                                             log_medincome2017 + log_population2017 + ptg_Unemployed2017 +
                                             ptg_white2017 + ptg_OwnerOccupied2017 + ptg_Mortgage2017 + ptg_Bachelor2017 + Gini2017 + 
                                             ptg_SavingsAccount2017 + FloodZone_Share_Dev_Area_tract):Post | 
            
            factor(visitor_home_tract_placekey) + factor(month) | 0 | visitor_home_tract, data = FL_visitor_1804_1910)

a1 = felm(log(1+visitor_count_home_tract) ~ (FloodDepth_ft_tract + FloodDepth_ft_tract:FloodInsTakeup + 
                                               
                                             FloodDepth_ft_location + NAICS_2digit_name + dummy_brands +
                                             log_medincome2017 + log_population2017 + ptg_Unemployed2017 +
                                             ptg_white2017 + ptg_OwnerOccupied2017 + ptg_Mortgage2017 + ptg_Bachelor2017 + Gini2017 + 
                                             ptg_SavingsAccount2017 + FloodZone_Share_Dev_Area_tract):Post | 
            
            factor(visitor_home_tract_placekey) + factor(month) | 0 | visitor_home_tract, data = FL_visitor_1804_1910)

a2 = felm(log(1+visitor_count_home_tract) ~ (FloodDepth_ft_tract + FloodDepth_ft_tract:I(FloodInsTakeup>0.5) + 
                                               
                                             FloodDepth_ft_location + NAICS_2digit_name + dummy_brands +
                                             log_medincome2017 + log_population2017 + ptg_Unemployed2017 +
                                             ptg_white2017 + ptg_OwnerOccupied2017 + ptg_Mortgage2017 + ptg_Bachelor2017 + Gini2017 + 
                                             ptg_SavingsAccount2017 + FloodZone_Share_Dev_Area_tract):Post | 
            
            factor(visitor_home_tract_placekey) + factor(month) | 0 | visitor_home_tract, data = FL_visitor_1804_1910)

a3 = felm(log(1+visitor_count_home_tract) ~ (FloodDepth_ft_tract + FloodDepth_ft_tract:I(FloodInsTakeup>0.5) + 
                                               
                                               FloodDepth_ft_location + NAICS_2digit_name + dummy_brands +
                                               log_medincome2017 + log_population2017 + ptg_Unemployed2017 +
                                               ptg_white2017 + ptg_OwnerOccupied2017 + ptg_Mortgage2017 + ptg_Bachelor2017 + Gini2017 + 
                                               ptg_SavingsAccount2017 + FloodZone_Share_Dev_Area_tract):Post | 
            
            factor(visitor_home_tract_placekey) + factor(month) | 0 | visitor_home_tract, data = FL_visitor_1804_1910_matched)

a4 = felm(log(1+visitor_count_home_tract) ~ (FloodDepth_ft_tract + 
                                               
                                               FloodDepth_ft_location + NAICS_2digit_name + dummy_brands +
                                               log_medincome2017 + log_population2017 + ptg_Unemployed2017 +
                                               ptg_white2017 + ptg_OwnerOccupied2017 + ptg_Mortgage2017 + ptg_Bachelor2017 + Gini2017 + 
                                               ptg_SavingsAccount2017 + FloodZone_Share_Dev_Area_tract):Post | 
            
            factor(visitor_home_tract_placekey) + factor(month) | 
            (I(FloodDepth_ft_tract*FloodInsTakeup*Post) ~ I(FloodDepth_ft_tract*log(SCI_Flooded_Distant)*Post)) | 
            visitor_home_tract, data = FL_visitor_1804_1910)



# IV first stage regression ----------------------------------------------------
IV <- lm(FloodInsTakeup ~ log(SCI_Flooded_Distant) +
                    FloodZone_Share_Dev_Area_tract + log(medincome2017) + log(population2017) + ptg_Unemployed2017 +
                    ptg_white2017 + ptg_OwnerOccupied2017 + ptg_Mortgage2017 + ptg_Bachelor2017 + Gini2017 + 
                    ptg_SavingsAccount2017,
                  data = tract_data)

# F stats for IV first-stage ------
cragg_donald(
  ~FloodZone_Share_Dev_Area_tract + log(medincome2017) + log(population2017) + ptg_Unemployed2017 +
    ptg_white2017 + ptg_OwnerOccupied2017 + ptg_Mortgage2017 + ptg_Bachelor2017 + Gini2017 + 
    ptg_SavingsAccount2017,  #Controls
  ~FloodInsTakeup,           #Treatments
  ~log(SCI_Flooded_Distant),  #Instruments
  data =tract_data
)





stargazer(a0, a1, a2, a3, a4, IV,
          
          type = 'text',
          #         type = 'latex',
          omit = c("NAICS", "location", "2017", "brands", "Share"),
          omit.stat = c('ser'),
          omit.table.layout = 'n',
          
          float = F,
          column.sep.width = "2pt",
          notes.append = F,
          font.size="small"
          
)


